Unidad 3: Álgebra de Mapas y Geoprocesamiento

Manipulando datos espaciales

  • La tabla de atributos
  • Herramientas de cálculo
    • Funciones y variables
  • Trabajo con las geometrías
  • Conversión de geometrías
  • Conversión de modelos: raster <-> vector

Paquetes utilizados en esta lección

library(sf)
library(tidyverse)
library(tmap)
library(raster)
library(mapedit)
library(magrittr)
library(stars)
library(units)

Actividad 1. Clases, atributos y valores

  • Cargue las capas ocurrencias y Mask del fichero data/gpkg/tapir.gpkg
tapir_mask <- st_read("data/gpkg/tapir.gpkg", layer = "mask")
tapir_ocu <- st_read("data/gpkg/tapir.gpkg", layer = "ocurrencias")
  • Explore los datos y determine qué tipo de dato es cada campo.
glimpse(tapir_mask)
Rows: 1
Columns: 2
$ Name <chr> NA
$ geom <MULTIPOLYGON [°]> MULTIPOLYGON (((-81.47795 -...

Para el caso de mask, solo tiene un campo y es tipo character

glimpse(tapir_ocu)
Rows: 155
Columns: 7
$ Especie    <chr> "Tpinchaque", "Tpinchaque", "Tpinchaque", "Tpinchaque", "Tp…
$ Confirmado <lgl> TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE,…
$ timestamp  <dttm> 2008-08-18 09:48:42, 2009-09-29 19:37:12, 2009-02-12 19:47…
$ Longitud   <chr> "78º 21' 24.12\" W", "77º 34' 24.13\" W", "77º 32' 54.13\" …
$ Latitud    <chr> "0º 40' 38.24\" N", "2º 15' 8.24\" N", "2º 11' 38.26\" N", …
$ zona_UTM   <chr> "18N", "18N", "18N", "18N", "18N", "18N", "18N", "18N", "18…
$ geom       <POINT [°]> POINT (-77.6433 0.67729), POINT (-76.42663 2.25229), …

La capa de ocurrencias presenta 6 campos y en su mayoría son characters, también hay logical y datetime.

  • ¿Que sistema de coordenadas tienen las capas?

La capa de mask, tiene un epsg: 4326, mientras que la capa ocurrencias, tiene un epsg: 4326.

  • ¿Qué tipos de operaciones se pueden hacer con esos tipos de datos?

Con la columna tipo logical se puede realizar operaciones booleanas, somo sumar los verdaderos. Con las fechas probablemente una diferencia de fechas o conteo por un período de tiempo (e.g. por meses), una operación similar se podría hacer por Especie.

Actividad 2. Cálculos avanzados: Funciones personalizadas

Usando las capas del fichero “data/gpkg/muestreo.gpkg”:

muestreo <- st_read("data/gpkg/muestreo.gpkg")
  • Agregar una columna que contenga las coordenadas X y Y en el SRC UTM-17S.
muestreo$X <- st_coordinates(muestreo)[, "X"]
muestreo$Y <- st_coordinates(muestreo)[, "Y"]

knitr::kable(head(muestreo))
id ta_media Altitud Forest_P_2010 Forest_P_2000 geom X Y
0 11.6 3513.165 22.5325 18.28994 POINT (763860.6 9677462) 763860.6 9677462
1 18.8 1899.906 43.8950 25.13772 POINT (746943.6 9598100) 746943.6 9598100
2 18.8 1937.868 26.9850 9.45509 POINT (749704.4 9644103) 749704.4 9644103
3 16.9 2144.949 37.0575 24.38037 POINT (732393.7 9612334) 732393.7 9612334
4 11.8 3073.430 37.7225 38.11243 POINT (675170.9 9598945) 675170.9 9598945
5 15.2 2740.981 36.6025 20.14110 POINT (779156.3 9688753) 779156.3 9688753
  • Agregar una columna nueva que contenga los pesos ficticios entre 70 y 120 kg, asignados de manera aleatoria a cada registro
muestreo$Weight <- sample(70:120, size = nrow(muestreo), replace = T)

knitr::kable(head(muestreo))
id ta_media Altitud Forest_P_2010 Forest_P_2000 geom X Y Weight
0 11.6 3513.165 22.5325 18.28994 POINT (763860.6 9677462) 763860.6 9677462 74
1 18.8 1899.906 43.8950 25.13772 POINT (746943.6 9598100) 746943.6 9598100 96
2 18.8 1937.868 26.9850 9.45509 POINT (749704.4 9644103) 749704.4 9644103 109
3 16.9 2144.949 37.0575 24.38037 POINT (732393.7 9612334) 732393.7 9612334 99
4 11.8 3073.430 37.7225 38.11243 POINT (675170.9 9598945) 675170.9 9598945 104
5 15.2 2740.981 36.6025 20.14110 POINT (779156.3 9688753) 779156.3 9688753 113
  • Agregar una columna nueva que contenga la zona UTM a la que pertenecen
GetUtmZone = function(x){
  centroid = st_centroid(x)
  longitude = st_coordinates(centroid)[,1]
  latitude = st_coordinates(centroid)[,2]
  zone_number = floor(((longitude + 180) / 6) %% 60) + 1
  zone_letter = ifelse(latitude >= 0,
                    'N', #if true
                    'S') #else
  return(paste0(zone_number, zone_letter))
}

muestreo$Zone <- st_transform(muestreo, 4326) |> st_geometry() |> GetUtmZone()

knitr::kable(head(muestreo))
id ta_media Altitud Forest_P_2010 Forest_P_2000 geom X Y Weight Zone
0 11.6 3513.165 22.5325 18.28994 POINT (763860.6 9677462) 763860.6 9677462 74 17S
1 18.8 1899.906 43.8950 25.13772 POINT (746943.6 9598100) 746943.6 9598100 96 17S
2 18.8 1937.868 26.9850 9.45509 POINT (749704.4 9644103) 749704.4 9644103 109 17S
3 16.9 2144.949 37.0575 24.38037 POINT (732393.7 9612334) 732393.7 9612334 99 17S
4 11.8 3073.430 37.7225 38.11243 POINT (675170.9 9598945) 675170.9 9598945 104 17S
5 15.2 2740.981 36.6025 20.14110 POINT (779156.3 9688753) 779156.3 9688753 113 17S
  • Agregue una columna nueva que contenga un identificador para los puntos que pertenecen al hemisferio sur y que caigan dentro del área que cubre la Demarcación hidrográfica Santiago
DHS <- st_read(
  "data/gpkg/muestreo.gpkg", layer = "Demarcacion hidrogafica Santiago",
  quiet = T
)

muestreo$dhs_sur <- st_intersects(muestreo, DHS, sparse = F) & muestreo$Zone == "17S"

knitr::kable(head(muestreo))
id ta_media Altitud Forest_P_2010 Forest_P_2000 geom X Y Weight Zone dhs_sur
0 11.6 3513.165 22.5325 18.28994 POINT (763860.6 9677462) 763860.6 9677462 74 17S TRUE
1 18.8 1899.906 43.8950 25.13772 POINT (746943.6 9598100) 746943.6 9598100 96 17S TRUE
2 18.8 1937.868 26.9850 9.45509 POINT (749704.4 9644103) 749704.4 9644103 109 17S TRUE
3 16.9 2144.949 37.0575 24.38037 POINT (732393.7 9612334) 732393.7 9612334 99 17S TRUE
4 11.8 3073.430 37.7225 38.11243 POINT (675170.9 9598945) 675170.9 9598945 104 17S FALSE
5 15.2 2740.981 36.6025 20.14110 POINT (779156.3 9688753) 779156.3 9688753 113 17S TRUE

Actividad 3. Edición de capas de líneas y polígonos

  • Cargue la librería mapedit
library(mapedit)

  • Ejecute la función mapedit::drawFeatures()

  • Acérquese a la ciudad de Cuenca, específicamente al Parque Calderón (Parque central de Cuenca)

  • Trace puntos de los 2 o 3 edificios representativos que que pueda reconocerlos fácilmente alrededor del Parque.

  • Trace un rectángulo que cubra hasta tres cuadras alrededor del parque calderón.

  • Trace las vías hasta dos cuadras alrededor del Parque.

  • Guarde el resultado en una capa, con SRC EPSG:4326 (Latitud y Longitud).

parque_calderon <- drawFeatures()

Edición de puntos, líneas y polígonos

mapdrawing

knitr::kable(head(parque_calderon))
X_leaflet_id feature_type geom
1622 marker POINT (-79.00394 -2.898353)
1628 marker POINT (-79.0053 -2.897426)
1634 marker POINT (-79.00459 -2.898053)
2017 rectangle POLYGON ((-79.00741 -2.9017…
2058 polyline LINESTRING (-79.0069 -2.897…
2082 polyline LINESTRING (-79.00678 -2.89…

  • Agregue una columna con el nombre o referencia de las vías dibujadas.

  • Explore el resultado y discuta sobre las propiedades de esta capa.

Para esto visualizaremos el mapa con tmap y señalaremos las ids de las vías para asignar el nombre que le corresponde:

tmap_mode("view")
filter(parque_calderon, st_geometry_type(parque_calderon) == "LINESTRING") |> 
  qtm()

Asignación de nombres a las vías

parque_calderon %<>%
  mutate(
    vias_id = case_when(
      `X_leaflet_id` == 2358 ~ "Gran Colombia",
      `X_leaflet_id` == 2082 ~ "Simón Bolívar",
      `X_leaflet_id` == 2058 ~ "Mariscal Sucre",
      `X_leaflet_id` == 2154 ~ "Presidente Córdova",
      `X_leaflet_id` == 2174 ~ "Juan Jaramillo",
      `X_leaflet_id` == 2238 ~ "Padre Aguirre",
      `X_leaflet_id` == 2266 ~ "Benigno Malo",
      `X_leaflet_id` == 2292 ~ "Luis Cordero",
      `X_leaflet_id` == 2316 ~ "Presidente Borrero",
      T ~ "Otras entidades"
    )
  )

Es un tipo de multi-geometría, es decir tiene una mezcla de polígonos, puntos y líneas. A pesar de que R lo puede procesar, otros programas como Qgis pueden necesitar que estás geometrías estén separadas.

  • Separe las entidades por geometrías y guarde como capas separadas en un fichero *.gpkg.
filter(parque_calderon, st_geometry_type(geometry) == "LINESTRING") |>
  st_write("data/gpkg/parque_calderon.gpkg", layer = "pc_lines")

filter(parque_calderon, st_geometry_type(geometry) == "POINT") |>
  st_write("data/gpkg/parque_calderon.gpkg", layer = "pc_points")

filter(parque_calderon, st_geometry_type(geometry) == "POLYGON") |>
  st_write("data/gpkg/parque_calderon.gpkg", layer = "pc_polygons")

Actividad 4. Geometrías derivadas

  • Cargue la capa puntos disponible en el fichero quinta_balzay.gpkg
quinta_balzay <- st_read("data/gpkg/quinta_balzay.gpkg")
  • Explore las propiedades de la capa y haga las preguntas necesarias para entender los atributos.
  • Genere una o varias capas de según corresponda. Si es necesario convierta esas geometrías a geometrías más complejas.

Para evitar que las geometrías se agrupen incorrectamente, hay que aplicar el argumento do_union = F en summarise()

Formación de geometrías a polígonos

polygons_qb <- quinta_balzay |>
  filter(feature_type == "polygon") |>
  group_by(id) |>
  summarise(do_union = F) |> 
  st_cast("POLYGON")

Formación de geometrías a líneas

lines_qb <- quinta_balzay |>
  filter(feature_type == "polyline") |>
  group_by(id) |>
  summarise(do_union = F) |> 
  st_cast("LINESTRING")

Extracción de puntos

points_qb <- quinta_balzay |>
  filter(feature_type == "point")
qtm(polygons_qb, fill = "id") +
  qtm(lines_qb) + qtm(points_qb)
  • Explore los resultados y comente lo que descubrió.

Los polígonos hacen referencia a los aularios y a las oficinas cerca de hídrica, las líneas a las vías dentros del campus balzay y los puntos a zonas importantes de la casa antigua de la quinta balzay (e.g. laboratorios).

Actividad 5. Raster a vector

  • Cargue la capa de Temperatura_3.tif y convierta las clases a polígonos
  • A partir de la capa Temperatura_2.tif, genere isolíneas cada 5 grados de temperatura.
at_3 <- raster("data/tif/Temperatura_3.tif")
at_2 <- raster("data/tif/Temperatura_2.tif")

A continuación se recorta las capa de temperatura con la de muestreo, para ahorrar tiempo en el procesamiento de las capas.

  • Isolíneas cada 5 grados
at_2_iso <- at_2 |> crop(muestreo) |> 
  cut(breaks = seq(-5, 35, 5), include.lowest = TRUE, right = FALSE) |>
  raster::rasterToContour() |>
  st_as_sf()
  • Polígonos de las clases
at_3_poly <- at_3 |> crop(muestreo) |> st_as_stars() |> st_as_sf(merge = T)
qtm(at_2_iso, lines.col = "level")
qtm(at_3_poly, fill = "Temperatura_3")
qtm(at_3_poly, fill = "Temperatura_3") +
  qtm(at_2_iso, lines.col = "level")

Geoprocesamiento y Álgebra de mapas

  • Geoprocesamiento (Vector)
    • Análisis de solape
    • Análisis de proximidad
    • Análisis de agregación
  • Álgebra de mapas (Raster)
    • Operaciones locales
    • Operaciones zonales
    • Operaciones globales

Actividad 5. Raster a vector

Una laguna de interés de un estudio se ecuentra en la siguiente coordenada: POINT(-79.27317 -2.85857). Usando las capas del fichero geoprocesos.gpkg, encuentre las respuestas a las siguientes preguntas:

Creamos el punto espacial e importamos la capa de geoprocesos.gpkg

laguna_point <- "POINT(-79.27317 -2.85857)" |> st_as_sfc(crs = 4326) |>
  st_as_sf() |> st_transform(crs = 32717)

water_bodies <- st_read("data/gpkg/geoprocesos.gpkg", layer = "cuerpos de agua") |>
  st_transform(crs = 32717)
  • ¿En un radio de 6 Km desde el punto, cuántos cuerpos de agua naturales existen?
st_is_within_distance(
  x = laguna_point, y = water_bodies, 
  dist = units::set_units(6,"km"),
  sparse=FALSE
) %>% sum()
[1] 126

Existen 126 cuerpoos de agua al rededor de los 6km de la laguna de interés.

  • ¿Cuál es el área que cubren todos ellas?
lag_buff <- st_buffer(laguna_point, dist = units::set_units(6, "km"))

lag_buff |>
  st_area() |> units::set_units("km^2")
113.0457 [km^2]

Tiene un área total de 113.05 km\(^2\)

  • ¿Cuál es el rango de temperaturas que existen en la zona de influencia? (Use la capa de Puntos de muestreo disponible en muestreo.gpkg)
qtm(muestreo) + 
  qtm(lag_buff)
muestreo[lag_buff, ] |> pull(ta_media) |> mean()
[1] 8

Tiene una media de 8°C

  • ¿Cuánto espacio del radio de 6 km, no está cubierto por cuerpos de agua?
st_difference(
  lag_buff, st_union(water_bodies)
) %>% st_area() %>% 
  units::set_units("km^2")
106.433 [km^2]

No están cubiertos por agua 106.4 km\(^2\)

Actividad 6. Álgebra de mapas

En un proyecto nuevo ya sea en QGIS o R cargue la capa temperatura_2.tif. Usando las herramientas necesarias genere una capa como la que se muestra en la imagen.

tm_shape(at_2) +
  tm_raster(
    palette = "-RdYlBu", 
    midpoint = 14, 
    style = "cont"
  )
at_levs <- cut(
  at_2, breaks = c(-Inf, -1, 5, 15, 22, 25, 50)
)
tm_shape(at_levs) +
  tm_raster(palette = "viridis", midpoint = 3)

Actividad 7. Álgebra de mapas 2

**Usando los recortes de Landsat 8 (clip_RT_LC08_L1TP_010063*), calcule el NDVI para toda la zona de estudio. El NDVI se calcula mendiante la ecuación:**

\[ NDVI = \frac{NIR - Red}{NIR + Red} \]

Importamos los datos y necesitamos las bandas 5 y 4

lc_f <- list.files(
  "data/tif", pattern = "clip_RT", 
  full.names = T
)

Cálculo de NDVI

ndvi_fun <- function(nir, red) {
  (nir - red) / (nir + red)
}

ndvi <- ndvi_fun(
  nir = raster(lc_f[4]), red = raster(lc_f[3])
)